function [S, SubKernel] = testrank( Kernel, Source, Target)
% Source / Target takes the form
% [ x_0 x_1 ]
% [ y_0 y_1 ]  from (x_0,y_0) to (x_1, y_1)


ntotal = length(Kernel);
nx = round(sqrt(ntotal));
ny = nx;

indSource = [];
indTarget = [];

nSourceCol = Source(1,2)-Source(1,1)+1;
if( Source(2,1) > 0 && Source(2,2) <= ny )
  indY = [Source(2,1) : Source(2,2)]';
elseif( Source(2,1) < 0 && Source(2,2) <= ny )
  indY = [(Source(2,1)+nx):ny,  1:Source(2,2)]';
elseif( Source(2,1) > 0 && Source(2,2) > ny )
  indY = [Source(2,1):ny,  1:(Source(2,2)-ny)]';
else
  error('Source dimension problem');
end
for ix = Source(1,1) : Source(1,2)
  indx = mod(ix-1, nx)+1;
  indSource = [indSource; indY + (indx-1) * ny];
end


nTargetCol = Target(1,2)-Target(1,1)+1;
if( Target(2,1) > 0 && Target(2,2) <= ny )
  indY = [Target(2,1) : Target(2,2)]';
elseif( Target(2,1) < 0 && Target(2,2) <= ny )
  indY = [(Target(2,1)+nx):ny,  1:Target(2,2)]';
elseif( Target(2,1) > 0 && Target(2,2) > ny )
  indY = [Target(2,1):ny,  1:(Target(2,2)-ny)]';
else
  error('Target dimension problem');
end
for ix = Target(1,1) : Target(1,2)
  indx = mod(ix-1, nx)+1;
  indTarget = [indTarget; indY + (indx-1) * ny];
end


S = svd( Kernel(indSource, indTarget) );
SubKernel = Kernel(indSource, indTarget);

if(0)
  indSource
  pause
  indTarget
  pause
  % imagesc(Kernel(indSource,indTarget))
  % pause
  % semilogy(S)
end
